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I. INTRODUCTION 

Considerable effort has been recently devoted to the 
understanding of magnetic thin films. The behaviour 
of magnetic insulators such as the transition-metal di- 
fluorides can be described in terms of short-ranged in- 
teraction models which makes the comparison with the- 
ory considerably simplerEJ. In particular, these materials 
can be epitaxially grown in very thin films to the point 
that specific theoretical concepts such as fitiite-size scal- 
ing close to a second order critical pointcrlj become, ex- 
perimentally verifiable. Such a study was pcrformedtJ for 
the (FeF2)n(ZnF2)m superlattice, where the magnetic in- 
teractions within a single FeF2 layer can be described in 
terms of a spin S — 2 bcc Ising model (with the different 
FeF2 layers sufficiently far apart that free boundary con- 
ditions can be assumed for each of them). The data for 
the thermal expansion coefficient a{T), which is exper- 
imentally observedo to be proportional to the magnetic 
contribution to the specific heat, show finite-size shifts of 
the critical point and rounding of the thermodynamic sin- 
gularity' in quantitative agreement with finite-size scaling 
theorytl. Subsequently, the thermal properties of (FeF2)„ 
(CoF2)n superlattices, where two different magnetic lay- 
ers interact, were studiedu. The thermal expansion co- 
efficient was studied as a function of temperature and 
of the layer thickness n. For n small, a{T) was found 
to show a single maximum, while for n larger, two ma:&- 
ima of a{T) as a funtion of temperature were observedu. 
Besides studying thermal properties, it is also possible 
to explore experimentally the magnetic properti^ of sin- 
gle monolayers through MoBbauer spectroscopyQ and to 
investigate the resulting order parameter profiles. 

In an attempt to provide a theoretical description of 
these layered magnetic systems beyond mean-field the- 
ory, we consider here as a simple toy model two coupled 
magnetic subsystems, where each subsystem contains n 
parallel layers of classical Ising spins, with 5=1 and 
S = \ respectively. We assume nearest neighbor cou- 
plings between the spins. For simplicity, we also assume 



that the coupling between spins within a layer is indepen- 
dent of S. We cannot expect with such a simple model 
to reproduce quantitatively any of the experiments men- 
tioned above, but we shall use our model as a means to 
obtain and to test scaling descriptions of the experimen- 
tally observed phenomena. Scaling should apply to more 
realistic situations. We shall thus write down the model 
in the form best suited for numerical treatment. 

Two main simplifications are employed. First, we 
work in two dimensions, considering layers of spin chains 
rather than the three-dimensional layers of films stud- 
ied experimentally. We expect that the scaling picture 
used to describe the critical behaviour can be applied 
in two as well as in three dimensions. Second, an ex- 
treme anisotropic limit is usedQ, where coupling constants 
between different layers are becoming very small while 
within a layer they become large. Then the task of cal- 
culating the thermal behaviour of the system amounts 
to studying the ground state properties of the quantum 
Hamiltonian 
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where a^--^ are the spin ^ Pauli jnatrices and S*^'^ are 
spin 1 matrices. It is well knownBfl that the critical be- 
haviour of this quantum chain is in the same universality 
class as the two-dimensional model of classical Ising spins 
described above (experimentiilly, this correspondence has 
recently been demonstrated^ for the dipolar-coupled 3D 
Ising ferromagnet LiIIoF4), but the numerical treatment 
of H is considerably easier than the corresponding cal- 
culation in the classical spin model using the transfer 
matrix. The critical behaviour of several coupled Ising 
systems witlLspin i in all subsystems had been investi- 
gated earlierBlHl. 

Let us explain the terms arising in H by making the 
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analogy with the two-dimensional model of classical Ising 
spins. The terms a^a^^-^ describe the interactions be- 
tween spins in different layers and the terms af describe 
the interactions within a single spin layer (and similarly 
for the S^'^). The coupling t plays the role of a tempera- 
ture. The S'-independence of the transverse field t reflects 
our assumption that the spin-spin coupling within a layer 
is spin-independent. Finally, 7 describes the coupling be- 
tween the two subsystems. The spatial coordinate £ cor- 
responds to the direction perpendicular to the magnetic 
layers. 

The following symmetries of H are immediate. First, 
H is invariant under the global spin reversal ct^ — > 
-c^^ 5^ ^ -5^ cr^ ^ cr^, S'' ~^ S"" . Those states which 
are invariant under this transformation are said to be 
even, all other states are said to be odd. Thus H is block- 
diagonalized into an even and an odd sector. Second, the 
spectrum of H is independent of the sign of 7, because 

(7) is changed into 7) through the similarity trans- 
formation tr^ — CT^, ^ S^. 

For each subsystem alone, that is for ^ ^ Q and n 00, 
there is a critical point at t = s withli3ll3 



= 1 



tc.i/n = 1.32587(1) 



(2) 



respectively. Varying k thus allows to change the ra- 
tio between the critical points tc,s in the S* = ^,1 sys- 
tems. Finally, C is a normalization constant which will 
be needed below in connection with the conformal invari- 
ance description of the spectrum of H at criticality. 

We are interested in the following observables,|Which 
will be studied through their quantum analogues .slJ The 
free energy of the two-dimensional classical spin model 
corresponds to the ground state energy Eq of H. Sim- 
ilarly, thermal averages < X > correspond to ground 
state expectation values {0|X|0). We also need the char- 
acteristic lengths ^1.2 of the spin-spin and energy-energy 
correlations (for r — )■ 00) 



G^{r) ^< cr(r)CT(0) > - < cr(r) >< a(0) > - e"''/*! 
G,{r) e(r)e(0) > - < e(r) >< e(0) > - e 
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where r is parallel to the individual layers. One has i^f 2 — 
-£'1,2 — Eq where Ei^2 are the energies of the first excited 
states in the odd and the even sector, respectively. 

The numerical technique used is completely standard, 
see Refs. for details. We use the Lanczos algorithm 
to find the first few lowest eigenvalues of H and the cor- 
responding eigenvectors. Finite-size scaling is then used 
to obtain estimates for the critical quantities which are 
then numerically extrapolated for n — > 00. 

This paper is organized as follows. In section 2, we 
discuss the phase diagram and comment on a subtlety in 
finite-size scaling. Section 3 describes the calculation of 
the surface critical exponents in 2D through conformal 
invariance techniques. In section 4, we present our results 
for the order parameter profiles. Finally, we give our 
conclusion in section 5. 



II. THE PHASE DIAGRAMS 

Our starting point is the experimental observation^ 
that the specific heat C as a function of the temperature 
will show one or two maxima depending on the thick- 
ness n. We therefore begin, with a consideration of this 
quantity. However, the explicit calculation of the second 
derivative —td'^Eo{t)/dt'^ of the free energy is cumber- 
some. To avoid this, recall the fluctuation-dissipation 
relation C ~ J2r^<i(^) together with the scaling form, 
which should be valid near criticality 



a(r) = r-2-^5(r/6) 



(4) 



where Xe = {1 — <y)lv and a,v are conventional critical 
exponents. Then, up to nonsingular background terms, 
the relation 



(5) 



should holdlla. Since we are here only interested in the 
leading critical behaviour, it is sufficient for us to con- 
sider the second gap £,2^ — E2 — Eq of H. The scaling 
behaviour of ^2 is simply related to the scaling of the spe- 
cific heat C and moreover the temperature dependence 
not too far away from the critical region of both ^2 (t) and 
C{t) should be qualitatively similar. Finalbi,j^2 is readily 
calculated through the Lanczos algorithmtflLl. We point 
out that the spin correlation length ^1 does not enter 
into the scaling form (||), because it couples to quantities 
which are odd under spin reversal while is evenllj. 

In figure [l] we show In ^2 as a function of t for different 
layer thicknesses n. 
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FIG. 1. The energy correlation length ^2 as a function of t 
and different layer thicknesses n for k = 4 et 7 = 1. The size 
of the layers is n = 2 ... 7 from bottom to top. 

We observe that for a very thin layer (n — 2, 3), there 
is only a single maximum present while two maxima de- 
velop for larger values of n. Comparing the location of 
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the maxima for n finite with the known values from (g) 
for their n — > oo limit, we see that the shift in the effective 
critical temperatures are quite large. Both maxima ap- 
pear to show a systematic build-up normally considered 
typical of a_thermodynamic singularity rounded by finite- 
size effectso. These observations, of one or two maxima 
depending on the value of n, large finite-size shifts of the 
pseudo critical temperatures and a rounding of the ther- 
modynamic-siagularity, are in qualitative agreement with 
experimeniu'liZl. 

Before we can make this conclusion however, one 
should realize that the models usually considered in the- 
oretical calculations and the superlattices studied exper- 
imentally are different. We shall refer to these as case A 
and case B, respectively. These cases differ in the way 
one goes from the finite system to an infinite one and it 
is only for the infinite system where a true phase tran- 
sition can occur. Consequently, the phase diagrams for 
cases A and B are different. (We reemphasize that in the 
following discussion we refer to two-dimensional systems 
while experiments are carried out in three dimensions.) 

A) We take a layer of n spins | and a layer of n spins 1. 
Periodic boundary conditions as explicitly written 
in are used. The phase diagram which results 
when n — » oo is given in figure 

B) In experimentsBSS, a procedure analogous to the 
following is used. One takes n spins ^ and n spins 
1 and repeats this double layer m times. Typically, 
n is small and fixed but m is large and formally, 
one should take a limit m ^ oo. This leads to the 
phase diagram in figure 
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(a) (b) 
FIG. 2. Phase diagrams for the two thermodynamic lim- 
its, for two-dimensional systems, (a) Two layers with n spins 
each, periodic boundary conditions and n ^ oo. (b) Super- 
lattice of m bilayers of n spins of each kind and m ^ oo. 

We first consider case A. Then, each of the two sub- 
systems can develop long-range order by itself. Conse- 
quently, the phase diagram (figure ^) will show four dif- 
ferent phases. There is a paramagnetic phase where the 
whole system is disordered, two distinct phases where 
either the spin i variables (< a >^ 0) or the spin 1 
variables (< S >^ 0) are ordered while the other subsys- 
tem is disordered and a ferromagnetic phase where the 



system is fully ordered. The transition lines are given 
by eq. (^) (full line for i and dashed line for ic.i)- In 
this case, the maxima observed in figure ^ should be in- 
terpreted as true thermodynamic singularities rounded 
by finite-size effects. Since we shall below concentrate on 
the properties of the order parameter close to the subsys- 
tem boundaries, we label these transitions by its surface 
critical propjerties, following the theory of surface phase 
transitional3. For the transitions from the paramagnetic 
phase to one of the partially ordered phases, one of the 
subsystems is still disordered and the order parameter of 
that subsystem which undergoes ordering will vanish at 
the boundary between the two subsystems. Along this 
line we have an ordinary transition. On the other hand, 
for the transitions from the partially ordered phases to 
the ferromagnetic phase one subsystem is already ordered 
which fixes the order parameter of the other subsystem 
at the subsystem boundary. Here wc have an extraordi- 
nary transition. At the meeting [Ppint of the transition 
lines there is a special transitiorJla. The scaling of the 
order parameter close to the subsystem boundaries is de- 
scribed by a different exponent than for the bulk, see 
Ref. ^ These local critical exponents are in 2D read- 
ily calculated using conformal invariance techniques, see 
section 3. 

For case B, corresponding to figure ^ however, the 
situation is different. Since each of the subsystems only 
contains a finite number of layers, the superlattice can 
only order as a whole. Thus the phase diagram contains 
a paramagnetic phase and a single ordered ferromagnetic 
phase. If t is sufficiently small, a layer of n spins may 
act as a giant spin and produce a strong thermal signal 
leading to a finite maximum of the specific heat or of 
related quantities. Since n is finite, however, there is no 
long-range order and the magnetic moments of each layer 
are independent of each other. Then the specific heat as 
a function of temperature will show two peaks, but only 
the one at lower temperatures will then correspond to a 
(shifted and rounded) phase transition and will develop 
a true singularity as m — > oo. Working in the framework 
of case B, it is misleading to call the location of the larger 
temperature maximum a (pseudo) critical point. 

From now on we always consider case A. Then, both 
maxima in ^2 can be interpreted as signalling a tran- 
sition. Also, we shall perform the subsequent scaling 
analyses just for the two-dimensional system, since the 
changes which might be needed in three dimcnsiotts, are 
immediate and discussed in detail in the literaturcBfi. 

How can one find the critical points from the finite- 
lattice data, when the Hamiltonians are more compli- 
cated and precise information on their location such as 
is not available a priori ? Practically, the transition 
points are located using phenomenoloKi|C|al renormaliza- 
tion as derived from finite-size scalingaB. Consider the 
quantity R{t;n) = n/^{t;n). Then finite-size estimates 
for the critical point tc can be found by solving for t the 
equation 
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R{t;n) = R{t;n + 1) (6) 

if n sufficiently large. A final value for tc is then obtained 
by extrapolating the resulting sequence for n —^ oo. Car- 
rying out this procedure, a further subtlety is encoun- 
tered as iUustrated in figure |[ 
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FIG. 3. Scaled inverse correlation lengths n/5i,2(t) as a 
function oft, for k = 4 and 7 = 1 and several layer thicknesses 
n. The critical points are labelled t° and corresponding to 
the ordinary and extraordinary transitions. The lower (up- 
per) curves correspond to (^2). 

Considering the finite-size scaling of the spin corre- 
lation length ^1 (lower curve), we see that the curves 
Ri{t;n) intersect close to the critical point t°. This is 
the cotLventional behaviour found e.g. in simple Ising 
modelsa'tl. For smaller values of i, vanishes exponen- 
tially fast with n, which refiects the ordering of the S = 1 
subsystem in this case. Thus in order to find tj;, the sec- 
ond gap $^2^ is the natural quantity to look at and in 
fact represents in the partially ordered phases the lowest 
physical excitation, just as does in the paramagnetic 
phase. Nevertheless, the curves R2{t;n) go through a 
minimum close to tc and will eventually touch each other 
in the n ^ 00 limit, but do not intersect. Although this 
is not in contradiction with the theory of finite-size scal- 
ing (at t = tc, (^) is strictly valid for n ^ 00 only), it is 
remarkable that at this point the conventional finite-size 
techniques are no longer applicable. In order to get an 
estimate of tc, one has to rely on locating a minimum of 
i?2 {t) or some other criterion. We stress that at t = t^ is 
the only phase transition occuring in the model for case 
B. 

This type of behaviour should be generic and although 
the example given does suggest that finite-size techniques 
may be fruitfully employed in analysing experimental 
data, it also shows that some care may be required. We 
shall see in the next section that in spite of the slightly 
unusual finite-size scaling, the spectrum of H at all these 



critical points is in full agreement with the conformal 
invariance predictions. 

To illustrate to what extent quantitative information 
about the critical behaviour can be extracted from our 
still relatively small systems {n < 7), we consider the 
determination of the correlation length exponent v. We 
look at the local maxima t^ax (n) Qir£2 near to the point 
t = t°. Finite-size scaling predictsoH that the tempera- 
ture shift 

I^t{n) =tl-tmax{n) ^ rT^/" , n ^ 00 (7) 
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FIG. 4. Determination of f from the finite-size scaling of 
tc{n) near the ordinary transition for n — 4. 

In figure ^, we show a log-log plot of At{n) vs. n 
and find that the asymptotic behaviour (|^) is already 
realised for small n in our toy model, although the tc{n) 
are not at all close to the n ^ 00 value t° , see figure 
From the slope in figure ^, we read off 1/ ~ 1.03, in good 
agreement with the exact result v — 1. We point out that 
also (SD) data for the thermal expansion coefficient for 
FeFe./ZnF2 superlattices were successfully analysed this 
wa5iD, although also in this case At(n) is largeEH, leading 
to V — 0.64(4) in agreement with the theoretical value 
V ~ 0.63 for the 3D Ising model. 



III. CRITICAL EXPONENTS AND CONFORMAL 
INVARIANCE 

We now describe the calculation of the critical expo- 
nents. Since we work with a two-dimensional classical 
spin model universalitjyi-elass, we can use conformal in- 
variance techniquesQcj l3 for that purpose. Here the use 
of the quantum Hamiltonian rather than the classical 
spin model becomes advantageous, since the surface crit- 
ical exponents which describe the local scaling close to 
the boundary between the two subsystems and in which 
we are mainly interested are easily obtained from the 



4 



low-lying excitation spectrum of H. We can thus avoid 
the cumbersome procedure of calculating first an average 
< X > and then subtracting from it its bulk contribution 
X\j in order to get the surface term Xg =< X > —Xb and 
then analyse its scaling behaviour. In the next two sub- 
sections, we shall first briefly collect the necessary back- 
ground knowledge and shall then apply it to the problem 
at hand. 



A. Ordinary and extraordinary transitions 

In two dimensions (and consequently also for quan- 
tum chains) conformal invariance specifies completely the 
scaling dimensions of all local observables. For a given 
model the first few exponents are very easily identified 
from the spectrum of the Hamiltpaian H. For free or 
fixed boundary conditions one haaH3 

= E,-Eo = n-\x, (8) 

where the exponents the local critical exponents 

which describe scaling near the boundary between the 
two subsystems and the index i labels the various scaling 
fields which occur in the model (usually, i = I corre- 
sponds to the order parameter and i = 2 to the energy 
density and higher gaps correspond to the scaling fields 
which generate correction-to-scaling terms) . The scaling 
of the gaps (||) goes with n and not with L because only 
half of the system is critical at either the ordinary or 
extraordinary transitions. However, in order to be able 
to apply (||) to the spectrum of a quantum chain such 
as (|l|), the normalization of H must be chosen such that 
energies and momenta are measured in the same units. 
One way of doing this is to recall that the surface scaling 
dimension of the energy density 

x,,s = 2 (9) 

which fixes the normalization C. Furthermore, once the 
normalization is fixed accordingly, the conformal alge- 
bra acts as a dynamical symmetry which determines the 
spectrum of H at criticality, viz. 

E, -Eo = -Lo + o(n-i) (10) 
n 

with Lq being one of the generators of the Virasoro alge- 
bra 

[Lj,Lk] = (j - k)Lj+k + ^(/ - j)5j+kfi (11) 

The universality class is determined byrtber-s^alue of the 
central charge c, for the 2D Ising modell3'c30 c = ^. 

Furthermore, the spectrum of H at criticality can be 
found from the representations of the Virasoro algebra, 
see Refs. ^,^0 for details. These representations are 
built from a highest weight state |A) which is defined 
through 



Lo|A)=A|A) , i,|A)=0 ifi>0 (12) 

and acts as the ground state for a certain representa- 
tion. Excited states are generated by acting with the 
L-j {j > 0) on |A). Now, the principle of unitarity 
of the underlying field theory restricts through the Kac 
formula the possible values of c and for each value of c 
only permits a finite number of possible values of A. For 
c — ^, the only possible values are A = 0,j^,^. This 
leads to the three unitary irreducible representations (0), 
(1/16) and (1/2) of the c — ^ Virasoro algebra. Now, 
the spectrum of H for thej-acdjnary and the extraordi- 
nary transitions is given byEjll3 

= (0) + (1/2) ordinary 

H^"'' = 2(0) extraordinary (13) 

where the trivial prefactor ir/n is suppressed. The factor 
2 for the spectrum at the extraordinary transition means 
that each level has the double degeneracy of the represen- 
tation (0). Combining these predictions with the formula 
(^) for the energy gaps, the critical exponents Xi can be 
read off. 

These predictioas— which had already been checked for 
the spin i beforeEj^ci], are fully reproduced in our model, 
in agreement with the expected universality. As an exam- 
ple, we take k = A and 7 = 1, but the results for the expo- 
nents do not depend on these parameters. The values for 
t at the ordinary and extraordinary transition are from 
eq. (1) t° = 5.30348(4) and = 1. The energy gap which 
is related through (||) to the exponent x^^s is the lowest 
gap in the even sector. Lattices with up to n = 7 were 
used. After extrapolatioiJj'Ea, we find C'°^ = 2.319(6) 
and C^''^ = 0.996(5) for the Cfd inary and the extraordi- 
nary transitions, respectiveljc^. In table |, we give the 
extrapolated estimates for the first four rescaled gaps 
Xi — {Ei ~ E[))n/{C,n) for the ordinary and extraordi- 
nary transitions together with the predictions following 
from conformal invariance. For the extraordinary transi- 
tion, all levels were found to be doubly degenerate in the 
n ^ oo limit but with an exponentially small splitting 
between pairs of levels for n finite. As should be expected 
from the algebraic construction of the spectrum of H, the 
differences between values of the Xi which belong to the 
same representation are integers. This reconfirms our de- 
termination of For the ordinary transition, we see that 
Xe,s = = 2 and xi — Xa-,s — Pi/v — 1/2 is the surface 
magnetization exponent, where j3i describes the scaling 
of the order parameter at the surface, mi ~ — T)^^ 
and V is the bulk correlation length exponentta. 

For the extraordinary transition, a little care is nec- 
essary in identifying the surface exponents. The order 
parameter is odd under spin reversal and the most local 
of all scaling operators. When this operator acts on the 
ground state, it creates the state with the lowest gap in 
the spectrum. Thus in our model x^n}s = = 0. On 
the other hand, in the literature (e.g. Ref. |l^ and refs. 
therein) , the extraordinary transition is defined with re- 
spect to those degrees of freedom which become critical in 
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the presence of a boundary which is already ordered. To 
read off the corresponding exponents, one should discard 
the double degeneracy of the spectrum, which is merely 
due to the ordering of the other subsystem. We then 

(2) 

have Xt^s = Xa^s = xi = 2_Xhis is in agreement with the 
expected scaling relationtZlES Xa^s — Pf' /v = 2 — a = 2. 



B. Special transition 

At the special transition, both subsystems become 
critical simultaneously. In addition, the Ising quantum 
chain has the peculiarity that the boundary coupling 7 
is marginal. The critical behaviour of the model can be 
described using previous results for the scaling behaviour 
of an Ising model with (semi-)infinite defect i±Hes,_sj£hich 
has been extensively studied for a long timeaiiil'Ej e3, see 
Rcf. ^ for a review. The local critical exponents depend 
continuously on the coupling 7. The mapping of coupled 
Ising layers to a 2D Ising model with a star-like config- 
uration of semi-infinite defect lifter was exactly derived 



for coupled spin ^ Ising modelsLl lj. The surface critical 



exponents can be read off the energy spectrumc3 

= E., -Eo = L-^2t:x,{-/) 



(14) 



provided that the normalization is fixed such that con- 
formal invariance is applicable. We find C from the re- 
quirement Xg^s — 1, which is also a necessary condition 
for the marginality of the coupling 7. 

The conformal theory is in this case more complicated 
than for the ordinary or extraordinary transitions. For 
the spin i Ising model, one can construct the Hamil- 
tonian spectrum either through j^on-unitary Virasoro 
generatora^j, Kac-Moody algebrasEi-er alternatively rely 
on boundary conformal field thcoryl23. Here we shall re- 
strict ourselves to a simple way to characterize the spec- 
trum. 

Taking the spin ^ case as a guide, we expect that for 
n large, the low-lying excitation spectrum QL.fi, can be 
recovered from the free fermion HamiltoniaiHI^il 



H 



r=0 



(15) 



where ni^"^ are fermionic number operators and A de- 
pends on 7. Non-universal terms which do not enter into 
the gaps are already subtracted. In general, for states 
in the even sector (with an even number of occupied 
fermionic states) and in the odd sector, there will be dif- 
ferent values Ao(7) and Ai(7), respectively. Now the 
lowest levels of H can be easily written down in terms 
of Aq.i. For example, the lowest exponents in the odd 
sector are 



1) 



1/. .0 1 
^(Ai + l 



2^0 



and the lowest exponents in the even sector are 



Xe,s = 1 



2Ao 
2Ao 



(16) 



(17) 



All these exponents correspond to conformal highest 
weight states. In addition, conformal invariance implies 
that if the exponent Xi of a highest weight state occurs 
in the spectrum, also Xi -\- k with k — 1, 2, 3 ... is present, 
with a known degeneracy which only depends on k (and 
which is 1 for the lowest two levels). From equations 
( |l6| , p7| ), the values of A04 for a given 7 are fcHfli^n While 
these are known exactly for the spin ^ casd^oO, these 
have be determined numerically for the case at hand. 

We first fix the normalization constant C, from the con- 
dition x^^s = 1- This condition means that the scaled 
lowest gap Aeven = nAEeven = ^LAE^^f.^ in the even 
sector should be equal to tt^, see (p^). In table ||, we 
give Af,yen for several values of 7. We find that within 
our numerical accuracy, its value is independent of 7 (the 
apparent deviation seen for 7 = 0.5 is an artifact from 
the extrapolation of our short sequences and should dis- 
appear if larger lattices could be taken into account) and 
conclude that the normalization C, is independent of 7. 
That is only to be expected fropi [eaclier results for spin 
i Ising models with defect line£2l'Lnl^. The final value of 
C, is taken from the values of 7 = 0.75 . . . and 7 = 0.87 . . . 
where convergence is best and we obtain C, — 0.5964(2). 

The numerical estimation for the higher gaps is made 
difhcult by (a) the relatively short sequences available 
{n < 7) and (b) the fact that for n finite, level cross- 



ings between different sequences occur. In table III we 
give the extrapolated results for the critical exponents 
Xi = L/{2Tr(){Ei — Eq) for several values of 7. When no 
information is given, our sequences did not converge reli- 
ably. We now want to compare these with the spectrum 
following from (p^. First, we use ( p6|Jl7| ) to determine 
Aq.i, which are also given in table iHlT Depending on 
the value of 7, it turned out to be numerically prefer- 
able to fix first Aq from ( [l7| ) and than use this value and 
the estimate of Xu to find Ai or alternatively determine 
Ai from the difference Xa',s —Xa-^si which is independent 
of Aq. The values of Ao,i were then used to calculate 
the other exponents which are listed in table 



[II as 'ex- 



pected'. When no error is given in these columns, the 
expected value is exact. 

We see that in general the extrapolated estimates for 
the higher gaps agree with the conformal invariance pre- 
diction to within a few per cent. A particular problem 
arises for 7 — 0.5, where the converge for is particu- 
larly slow. In that case, we are not able to sensibly to 
specify accuracies for Ao,i and the correspondence be- 
tween the 'numerical' and the 'expected' data is more 
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qualitative. The situation here could only be improved 
by going to larger lattices. On the other hand, for the 
other values of 7, we obtain a nice agreement between 
the 'numerical' data and the 'expected' free fermion spec- 
trum. We point out that the beginning of several confor- 
mal towers (that is, with Xi also Xi + 1 and even Xi + 2 are 
found in the spectrum) is observed. The fact that this 
level spacing comes out correctly is a further confirma- 
tion of our determination of the normalization constant 
(. On the other hand, we have not been able to go suf- 
ficiently high in the spectrum to check the degeneracies 
of the excited states. 

We see that the scaling behaviour of our model is de- 
scribed in terms of a free fermion system. This should be 
expected on the basis of universality, although this free 
fermion description of a spin 1 model is not at all obvi- 
ous from the lattice formulation. Nevertheless, there is 
an important distinction with respect to the spin ^ case. 
Recall that the value of 7 is the same at both subsys- 



tem boundaries 
would have foundEJ 
our model, see tabic 



e coupled two spin i systems, we 
Aq = 0, which is not the case in 



II] 



IV. ORDER PARAMETER PROFILES 

So far, we have calculated the critical exponents which 
describe the scaling of observables close to the subsys- 
tem boundary. We now ask for the form of the order 
parameter profiles close to that interface. 



A. Generalities 

The calculation of the order parameter on a finite lat- 
tice poses a conceptual problem. The natural candidate, 
< M >= (0|A/|0) = on any finite lattice. This diffi- 
culty can be overcome by first introducing a small mag- 
netic field h, calculating < M > in the presence of h, 
take the infinite system limit L — > 00 and only then let 
/i — > 0. In practice, rather than performing numerically 
this double limit, the following trick which goes back to 
Yang is used. In the ordered phase(s), the ground state 
is already on a finite lattice almost degenerate, where 
the energy splitting decreases exponentially with L. In- 
troducing an infinitesimal magnetic field h into H and 
working within degenerate first-order perturbation the.- 
ory in h, the order parameter on the site i is given bjtj 



m{i) = (l|M(i)|0) 



(18) 



where |0) and |1) are the lowest eigenstates in the even 
and odd sectors, respectively. For our model (|^), the 
magnetization operator M{i) is 



M{i) 



so that M{i) is normalised such that \M{i)\ < 1 for all 
sites. It is well knowriEj that the finite-lattice order pa- 
rameter calculated from (|l^) has the correct scaling be- 
haviour. The dependence of the order parameter profiles 
on S deep in the ordered phase has also been studied23. 

Practically, for the computation of the eigenvectors 
|0),|1), it is not necessary to store all the intermediate 
Lanczos vectors. This can be avoided by running the 
Lanczos algorithm twice, where the first pass furnishes 
the weights by which the intermediate vectors contribute 
to |0), |1) and in the second, pass the eigenvectors them- 
selves can be accumulated^. 



B. Ordinary and extraordinary transitions 

Before presenting our results for the order parameter 
profiles at the various transitions, leLjUa-adapt the predic- 
tions from finite-size scaling theoryou'EEl to the situation 
at hand. We are interested in the local order parameter 
m(i) rather than the full magnetization m = 
One should distinguish whether m{i) is measured far 
away or close to the subsystem boundary. In the first 
case, when the site i is well in the bulk, we expect 



m(i) 



'M 



2i-l 
4n 



(20) 



spm 2 region 
spin 1 region 



(19) 



In the secoud case, when i is close to the boundary, we 
should haveE3 

m(i) a-^- (^- j Mi — (21) 

Here, Xa- and the bulk and surface critical ex- 

ponents calculated in the previous section, M and M 
are scaling functions, i = 1, 2, . . . , 2n is measured from 
the left boundary of the 5=1 subsystem, n is the layer 
thickness and a the lattice constant. We point out that 
the arguments of the two scaling functions are different .EZl 
In the first case, the scaling is such that the total system 
size is kept fixed and the lattice constant a — > 0, while 
in the second case, a is kept fixed and the system size 
L = 2n ^ 00. 

These predictions are confirmed by our numerical re- 
sults. Consider first the ordinary transition. Again, we 
take K = 4 and 7 = 1 as an example. In figure ||, we 
rescale our magnetization profiles according to (^^ with 
Xa = 1/8 and we see that indeed for the portion of the 
lattice which is far enough for the subsystem boundaries, 
a data collapse occurs even for the small lattices consid- 
ered here. Also, we see that in the immediate vicinity 
of the subsystem boundaries, the scaling description ([2^) 
no longer applies. Similar plots for the bulk scaling can 
be obtained for the other transitions but will not be pre- 
sented here, but see figure || below. 
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1^ 

FIG. 5. Order parameter for the ordinary transition, 
scaled with the bulk exponent Xa ~ 1/8 for k = 4 and 7 = 1 
as a function of ^ = {2i — l)/4n. The regions of spin 5 = 1 
and S — ^ are indicated and the boundary between them is 
shown by the dotted line. The symbols correspond to n = 4 
(diamonds), n = 5 (triangles), n = 6 (squares) and n — 7 
(circles). 



a non-magnetic substrate show two-dimensional critical 
behaviour for layer thicknesses of less than about two 
monolayers and cross over to three-dimensional critical- 
ity for only slightly thicker layersH. 

A similar behaviour is also found for the extraordinary 
transition. However, as already mentioned in discussing 
the spectra, it is sensible to distinguish two 'order pa- 
rameters'. These are 



m^^^i) = (l|M(z)|0) 

m(2)(i) = (l|M(z)|0') ~ {l'\M{t)\0) 



(22) 



where |0') and |1') are the first excited states in the even 
and odd sectors and the approximate equality between 
the two forms for m^^-* holds up to terms exponentially 
small in L. Note that here the ordering at the subsystem 
boundary is provided through the subsystem already in 
its ordered phase for n oo and not through fixing the 
spins at the boundary. The profile for m^^\ where the 
surface exponent 2;^^], = 0, is shown in figure ^. Again, 
we see that for the S ~ ^ subsystem, we have a data 
collapse according to for the first two monolayers 
next to the boundary and for larger values of i, the is 
a rapid crossover toward the bulk scaling (pO|). Since 
the S — I subsystem is ordered, finite-size effects are 
exponentially small there. 



S=l/2 




S=l/2 



FIG. 6. Scaled profile of the order parameter for the or- 
dinary transition, for k = 4 and 7 = 1 as a function of 
jj, = i — n — ^. The inset shows the location of the tran- 
sition point in the phase diagram. The correspondence of the 
symbols to the layer thickness n is the same as in figure ^. 

In figure ^, we display the local scaling of the order 
parameter close to one of the subsystem boundaries ac- 
cording to (pl|). For the ordinary transition, x^.s = 1/2 
and we set a = 1. We see that in the first two monolayers 
around both sides of the interface, the data collapse onto 
the scaling form (pi]). However, going beyond the first 
two monolayers, it is apparent that there is a crossover 
towards the bulk scaling form (pO[). It is apparent that 
the surface scaling only occurs in a very thin layer close 
to the boundary. We remark that this is consistent with 
experimental observations that thin magnetic layers on 



0.8 




FIG. 7. Profile for the order parameter m'^' at the extraor- 
dinary transition for k = 4 and 7 = 1 as a function of fj,. The 
inset shows the location of the transition point in the phase 
diagram. The correspondence of the symbols to the layer 
thickness n is the same as in figure ^ 



C. Special transition 

This case is of particular interest, since the exponent 
Xct,s does depend on 7. It is therefore interesting to ask 
whether the profiles are affected by changing 7 as well. 
The bulk scaling behaviour ( ^0| ) with Xo- = 1/8 is recov- 
ered as in the other transitions. 
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FIG. 8. Scaled profiles for the order parameter at the spe- 
cial transition for several values of 7 as a function of [s,. 
The dotted curves with full symbols correspond to 7 = i, 
the fuU curves with open symbols to 7 = 0.877111 and the 
dash-dotted curves with open symbols corresponds to 7 = 2. 
The correspondence of the symbols to the layer thickness n is 
the same as in figure ^. 

In figure ^, we show for three values of 7 the local 
scaling of the order par ameter, where the values of Xa^s 
are taken from table III . On both sides of the boundary, 
we find a data collapse according to (|T]) for the first few 
boundary layers and for l arg er values of i, a crossover 
towards the bulk scaling (|20|). In addition, we see that 
the form of the scaling function M. does depend on 7. For 
7 ~ 0.87, we have Xa,s ~ = 1/8 (see table III) and the 
distinction between local and bulk scaling is somewhat 
washed out. 

Concerning the shape of the scaling function A^(/x), 
we see that for the special transition, it can be a non- 
monotonous function of pt. For the ordinary and the 
extraordinary transitions, however, it is a monotonous 
function of /I and this holds independently of the value 
of 7. For the special transition, A^(/I) is only monotonous 
if 7 ~ Kc- Qualitatively, this can be explained as follows. 

First, if 7 ~ Kc, the boundary coupling takes the ef- 
fective mean value which smoothly interpolates between 
the two different regimes. Since both systems become 
critical simultaneously at the special transition, the scal- 
ing functions simply interpolates smoothly between the 
values of the magnetization finite-size scaling amplitudes 
in the two subsystems. Since these amplitudes are differ- 
ent, even for Xa^s{l) = 1/8, the scahng function A^(pt) 
will not become a constant. Second, consider 7 > k^- 
Then the spins on both sides of the boundary are more 
strongly coupled together than two spins in either sub- 
system. Since the 5* = state in the spin 1 subsystem 
does not contribute to the energy, this leads to an en- 
hancement of states where the boundary spins on both 
sides are up. Indeed, we checked that already for 7 = 4, 



the local order parameter on both sides of the boundary 
is close to saturation. Thus, we have a large value of m(i) 
close to the boundary which then falls back to an average 
value for each of the subsystem, in agreement with fig- 
ure ||. Finally, for 7 < Kc, the spins on both sides of the 
boundary are more weakly coupled than average spins. 
This favors states with a smaller value of m(i) close to 
the subsystem boundary. 

On the other hand, for the ordinary or the extraor- 
dinary transition, one subsystem is much more ordered 
than the other one. If 7 is large, the first spin across the 
boundary is strongly aligned with the spins of the more 
ordered subsystem and if 7 is small, the coupling of the 
first spin to the more ordered subsystem is reduced. This 
leads to an effective translation of the order parameter 
profile without affecting its form. 



D. Magnetization profiles and conformal invariance 

In 2Z?, conformal invariance states that the profile of a 
local scaling operator ip with bulk scaling dimension x^ 
is on an infinitely long strip of finite width L and with 
the same type of boundary conditions on both sides ifend 
in particular for free boundary conditions) given byE3 



(23) 



where v measures the position across the strip and is 
a non- universal constant. The scaling function for mixed 
boundapz. conditions is also known for minimal conformal 
thcorieso. This result only depends on the transforma- 
tion properties of the scaling operator tp. Furthermore, 
this result carries over to the profiles on quantum chains. 

When we try to apply this to the order parameter at 
the ordinary transition, we should find Aa- = due to 
symmetry. However, the finite-lattice estimates for m(i) 
obtained from eq. (|l^) above involved an infinitesimal 
magnetic field h, which (a) invalidates the above symme- 
try argument and (b) leads to a new effective exponent 
x^p = Xa- — Xa.s- This is seen as follows. From our numer- 
ical data, we have found the scaling form ( pO| ) 

ruLiz) = L-""' M{n) , ^i = z/L (24) 

On the other hand, close to the boundary, the order 
parameter should show surface finite-size scaling ^ 
L"^"'". This implies for the scaling function, e.g. Ref. |l^ 



X(/.) ^ 







(25) 



from which be identified. Eq. ( |25| ) was alse, con- 

firmed within the framework of the e expansiorcll and 
for the Ising quantum chain with an apeijiadic modula- 
tion generated by the Fredholm sequence.eS For the 2D 
Ising model, we remark that also in the presence of a 
small surface magnetic fieldcj hi the spatial dependence 
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of the magnetization near to the surface scales as z^/^, 
which is the same as obtained from eq. (^). 

We now try to extend ( p5| ) , derived for pnall values of 
/i only, to larger values of /i. If we accepia that the es- 
timate eq. for the spontaneous magnetization trans- 
forms covariantly under conformal transformations, we 
can combine eqs. (^^5|). Then the exact finite-sjize scal- 
ing function at the ordinary transition would beCj 

M{p) = A„ {sm2-K^f (26) 

taking into account that for our model, only the section 
< < ^ is actually critical at the ordinary transition 

for K > Kc- 




0.4 I 1 1 1 1 1 0.95 I ' ' ' ' ' I ' ' ' ' 

0.1 0.2 0..1 0.4 0.5 0.5 0.6 U.7 U.8 U.9 I 0.5 0.6 0.7 0.8 0.9 

n 

FIG. 9. Comparison of the conformal invariance prediction 
for the scaled order parameter profiles m{i)n^^^ . Only the 
part of the system which is critical is shown. The graphs give 
for (a) the ordinary transition, (b) the extraordinary transi- 
tion with m'^' and (c) for m^-^^ at k = 4 and 7 = 1. The 
correspondence of the symbols to the layer thickness n is the 
same as in figure ^. 

In figure |9^, we compare ( p6| ) to the numerical data. 
First, we observe a data collapse from several system 
sizes onto a single curve. Second, the form of the scaling 
function agrees nicely with eq. ( P6[) . The same resul|t-|has 
also been found for the Hilhorst-van Leeuven modelciS. 
Since for the 2D Ising model, (|2^) remains unchanged, 
even in the presence of a small surface magnetic field, o 
Ai should be independent of a small hi for the ordinary 
transition. 

Let us compare the finite-size scaling functions for the 
profiles coming from (^) and ( |2^ ) . The first one is based 
on a continuum description of the profile in the half- 
infinite systerp-jwhich is then conformally transformed 
onto the stripES. Very close to the boundary, a contin- 
uum description may no longer be applicable. Indeed, 
for unitary conformal theories such as the Ising model, 
the critical exponents x^p > and the profile as it stands 
will diverge at the boundary, in disagreement with exist- 
ing numerical data. On the other hand, the second form 
is constructed to be consistent with both bulk finite-size 
scaling deep inside the system and surface finite-size scal- 
ing close to the boundary. In order to match this with 
the functional form required from conformal invariance. 



it is necessary to assume that the exponent Xa- —x^^s gov- 
erns the scaling of the matrix element (^8|) (calculated in 
an infinitesimal magnetic field which breaks global sym- 
metry) used to estimate the finite-size order parameter, 
rather than the conventional order parameter scaling di- 
mension > 0. While this approach certainly is in 
agreement with the numerical data for the whole strip, it 
is not yet understood how the above dimensional trans- 
mutation can be explained. 

In addition, we find that the same functional form also 
describes the order parameter profiles for the extraordi- 
nary transition, as shown in figure ^. The numerical data 
are again consistent with scaling (note that the overall 
scale in figure ^ is abount an order of magnitude larger 
than in the other two cases). Specifically, we find from a 
fit 

{A^(sin27r^)3/8 ^ A„ ~ 0.80 ord. 
A^^^(sin27r^)-i/i6 ^ ^(^i) 0.98 ex., m(i) 
Ai-^^(sin27r^)9/8 , A^^^ - 0.49 ex., m(2) 

(27) 

Tentatively, the profile for vnS^^ (which is not given by 
( p3| ) because the 'ordered' subsystem is still finite, see 
eq. (^2|)) can be explained as follows. This order param- 
eter is sensible to the ordering which occurs at the or- 
dinary transition. At the extraordinary transition, these 
degrees of freedom have become massive and thus have 
a short effective correlation length ^e//- Then the fluc- 
tuating spins would see a fixed boundary on one side 
but because f ^ L the other boundary should appear 
as open. Hpagver, for mixed boundary conditions, it is 
known thatEju Xa.s ~ jq- Then —x^p = Jq ^ = ^jg ^"'^ 
agreement with the numerical data. 

Finally, the above argument does not reproduce the 
profiles for the special transition. This is due to the 
fact that the, two-point correlation functions are more 
complicated^ than tla*. simple power-law form which un- 
derlies the derivationESi of (p3|). 

V. SUMMARY 

We have studied the transitions arising in a pair of 
magnetic layers, coupled through short-ranged interac- 
tions and described by Ising models. This study was 
motivated by ongoing experiments on similar systems. 
We have found the variation of the 'specific heat' with 
the temperature and studied the scaling of the order pa- 
rameter profiles at the phase transition points. Our aim 
was to check out a scaling analysis which should also be 
applicable to experimental data in 3D. 

We have reemphasized that the systems usually stud- 
ied in experiments and the models best suited for theo- 
retical analysis are not completely identical and care is 
needed in the comparison of the two, as exemplified in 
the two phase diagrams in figure 0. 
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To simplify the theoretical analysis, we used a two- 
dimensional|-™edel of layers of spin chains, although the 
experimentfO^ESfl involve two-dimensional layers stacked 
in the third dimension. This was for us no serious dis- 
advantage, since the scaling arguments used here jCaii be 
extended to three dimensions in a well-known waytTH. In 
addition, conformal invariance techniques could be used 
to simplify the calculation of the surface critical expo- 
nents we were interested in. This calculation reconfirmed 
the expected universality of the surface critical exponents 
found for the ordinary, extraordinary and special transi- 
tion (where in the 2D Ising model, the surface coupling 
is marginal). 

Our results on how the 'specific heat' depends on 
the layer thickeess n are qualitatively the same as seen 
experimentallycl. Our results also suggest that if the 
layer thickness can be precisely controlled, one might 
get a trade-off in no longer having to achieve a very 
fine temperature control and still being able to measure 
fluctuation-dominated critical exponents to good accu- 
racy. 

Our study for the order parameter profiles was moti- 
vated by the existing experimental techniques,to measure 
the magnetic moments of a single monolayeitl. While for 
an infinite system, the critical point magnetization van- 
ishes, for finite lattices a non-trivial finite-size scaling 
behaviour of the order parameter profiles is obtained. 
We found two types of finite-size scaling at the transi- 
tions. Far away from the subsystem boundaries, we get 
a bulk finite-size scaling (|2^) governed by the bulk expo- 
nent = P/v ~ 1/8. In the immediate vicinity (and 
on both sides) of the boundary, however, the order pa- 
rameter finite-size scaling ( pl| ) is governed by the surface 
critical exponent Xa^s = l^i/i^ whose value depends on 
the type of the surface transition. 

For the ordinary and the extraordinary transitions, the 
data for finite-size scaling of the order parameter profile 
scaling function M (measured in an infinitesimal mag- 
netic field) suggest in 2D an exact scaling function from 
consistency with bulk and surface finite-size scaling and 
with conformal invariance. It remains a challenge to de- 
rive a similar result for the special transition. It is not 
yet clear how to address the problem of calculating the 
surface scaling function M.. In any case, a continuum 
approach, which underlies the conformal invariance ar- 
guments used for the determination of M , does not seem 
to be feasible in that case. 

It is a pleasure to thank R. Camley, F. Igloi and L. Tur- 
ban for useful discussions. 
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1 


0.4994(6) 




1/2 




2 


2 


2 


1.496(7) 




3/2 




3.01(2) 


3 


3 


2 




2 




3.95(6) 


4 


4 


2.505(9) 




5/2 




4.9(2) 


5 




TABLE I. Conformal spectrum 


of the surface exponents Xi 


for the ordinary and extraordinary transitions at 


K = 4 and 


7 


= 1. The numbers in brackets give the estimated uncertainty 


in the last given digit. 







7 0.5 0.65 0.754222 0.877111 1 

Aever, 1.945(3) 1.889(7) 1.8731(8) 1.8745(5) 1.88(2) 

TABLE IL Scaled and extrapolated lowest gap A^vm = ttC in ths even sector at the special transition for several values of 
7. The numbers in brackets give the estimated extrapolation error in the last digit. 





7 = 0.5 


7 = 


0.754222 


7 = 


0.877111 


7 = 1 




i 


numerical 


expected 


numerical 


expected 


numerical 


expected 


numerical 


expected 


1 


0.231(1) 


0.20 


0.1436(6) 


0.144(2) 


0.1103(7) 


0.1103(5) 


0.0841(5) 


0.084(5) 


2 


0.971(3) 


0.91 


0.999(1) 


1 


1.000(1) 


1 


1.00(2) 


1 


3 


1.034(5) 


1 


1.10(3) 


1.072(4) 


1.113(6) 


1.1103(5) 


1.11(2) 


1.084(5) 


4 


1.095(3) 


1.20 


1.12(2) 


1.144(2) 


1.168(5) 


1.168(2) 


1.232(5) 


1.232(5) 


5 


1.72(3) 


1.70 


1.94(1) 


1.94(1) 


1.92(1) 


1.92(1) 


1.9(1) 


1.78(4) 


6 


1.92(3) 


1.91 


1.995(5) 


2 


2.00(2) 


2 


2.0(1) 


2 


7 


1.99(2) 


2 


2.03(3) 


2.06(1) 


2.11(1) 


2.08(1) 


2.1(1) 


2.084(5) 


8 




2.20 


2.11 


2.072(4) 




2.1103(5) 


2.18(3) 


2.22(4) 


9 




2.30 


2.145(8) 


2.144(2) 




2.168(2) 




2.232(5) 


Ao 


0.15 




0.030(5) 


0.040(5) 


0.11(3) 




Ai 


0.36 




0.464(2) 


0.529(1) 


0.574(4) 



TABLE III. Conformal spectrum of the scaling dimensions ^^(7) at the special point t = 1, k = 0.754222. The values of 
Ao,i used in comparing with the free fermion Hamiltonian (^) are also given. The numbers in brackets give the estimated 
uncertainty in the last given digit. 
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